#include <boost/random.hpp>
#include <kuai/misc/xyz.hpp>
#include "test.hpp"

namespace kuai {

	void TestXYZJob::doJob(std::ostream& out_log, std::ostream& err_log, StringMap& par) {
		size_t n = getParameter(par, "n", size_t(50000));
		RealNumber cutoff = getParameter(par, "cutoff", RealNumber(0.2));
		RealNumber box  = getParameter(par, "box", RealNumber(1));
		size_t seed = getParameter(par, "seed", size_t(0));
		if (seed == 0) {
			seed = size_t(::time(NULL));
		}
		

		boost::mt19937 gen(seed);
		boost::uniform_real<RealNumber> dist(0, box);
		boost::variate_generator<boost::mt19937&, boost::uniform_real<RealNumber> > die(gen, dist);
		std::vector<XYZ> data(n);
		for (std::vector<XYZ>::iterator 
			it = data.begin(); it != data.end(); ++it) 
		{
			it->x = die();
			it->y = die();
			it->z = die();
		}

		test1(data, cutoff, out_log);
		test2(data, cutoff, out_log);
		test3(data, cutoff, out_log);
		test4(data, cutoff, out_log);
		test5(data, cutoff, out_log);
		test6(data, cutoff, out_log);
		test7(data, cutoff, out_log);
		test8(data, cutoff, out_log);


		// testC1(data, cutoff, out_log);
		testC2(data, cutoff, out_log);
		testC3(data, cutoff, out_log);
		testC4(data, cutoff, out_log);
		testC5(data, cutoff, out_log);
		testC6(data, cutoff, out_log);
	}

	void TestXYZJob::test0(const std::vector<XYZ>& data, RealNumber cutoff, std::ostream& out_log) {
		Time t1 = now();
		size_t hits = 0;
		for (std::vector<XYZ>::const_iterator 
			i = data.begin(); i != data.end(); ++i)
		{
			for (std::vector<XYZ>::const_iterator 
				j = i + 1; j != data.end(); ++j)
			{
				XYZ delta(*i-*j);
				if (delta.abs() < cutoff) {
					hits += 1;
				}
			}
		}

		Time t2 = now();

		out_log << "Test 0: Standard, Use iterator and abs\n"
				"hits " << hits << " and it costed " << t2-t1 << std::endl;
	}

	void TestXYZJob::test1(const std::vector<XYZ>& data, RealNumber cutoff, std::ostream& out_log) {
		Time t1 = now();
		size_t hits = 0;
		for (std::vector<XYZ>::const_iterator 
			i = data.begin(); i != data.end(); ++i)
		{
			for (std::vector<XYZ>::const_iterator 
				j = i + 1; j != data.end(); ++j)
			{
				XYZ delta(*i-*j);
				if (shorter(delta, cutoff)) {
					hits += 1;
				}
			}
		}

		Time t2 = now();

		out_log << "Test 1: Use iterator and shorter\n"
				"hits " << hits << " and it costed " << t2-t1 << std::endl;
	}

	void TestXYZJob::test2(const std::vector<XYZ>& data, RealNumber cutoff, std::ostream& out_log) {
		// 
		Time t1 = now();
		size_t hits = 0;
		size_t n = data.size();
		for (size_t i = 0; i < n; ++i) {
			for (size_t j = i+1; j < n; ++j) {
				XYZ delta(data[i]-data[j]);
				if (shorter(delta, cutoff)) {
					hits += 1;
				}
			}
		}

		Time t2 = now();

		out_log << "Test 2: Use index and shorter with buffer data.size()\n"
				"hits " << hits << " and it costed " << t2-t1 << std::endl;
	}

	void TestXYZJob::test3(const std::vector<XYZ>& data, RealNumber cutoff, std::ostream& out_log) {
		Time t1 = now();
		size_t hits = 0;
		for (size_t i = 0; i < data.size(); ++i) {
			for (size_t j = i+1; j < data.size(); ++j) {
				XYZ delta(data[i]-data[j]);
				if (shorter(delta, cutoff)) {
					hits += 1;
				}
			}
		}

		Time t2 = now();

		out_log << "Test 3: Use index and shorter without buffer data.size()\n"
				"hits " << hits << " and it costed " << t2-t1 << std::endl;
	
	}

	void TestXYZJob::test4(const std::vector<XYZ>& data, RealNumber cutoff, std::ostream& out_log) {
		Time t1 = now();
		size_t hits = 0;
		for (size_t i = 0; i < data.size(); ++i) {
			for (size_t j = i+1; j < data.size(); ++j) {
				XYZ delta(data[i].x-data[j].x, data[i].y-data[j].y, data[i].z-data[j].z);
				if (shorter(delta, cutoff)) {
					hits += 1;
				}
			}
		}

		Time t2 = now();

		out_log << "Test 4: Standard with explicit ctor\n"
				"hits " << hits << " and it costed " << t2-t1 << std::endl;
	}

	void TestXYZJob::test5(const std::vector<XYZ>& data, RealNumber cutoff, std::ostream& out_log) {
		Time t1 = now();
		size_t hits = 0;
		for (std::vector<XYZ>::const_iterator
			i = data.begin(); i != data.end(); ++i) {
			for (std::vector<XYZ>::const_iterator
				j = i+1; j != data.end(); ++j) {
				XYZ delta(*i, *j);
				if (shorter(delta, cutoff)) {
					hits += 1;
				}
			}
		}

		Time t2 = now();

		out_log << "Test 5: Same to 4 with iterator\n"
				"hits " << hits << " and it costed " << t2-t1 << std::endl;
	}

	void TestXYZJob::test6(const std::vector<XYZ>& data, RealNumber cutoff, std::ostream& out_log) {
		Time t1 = now();
		size_t hits = 0;
		const size_t n = data.size();
		for (size_t i = 0; i < n; ++i) {
			for (size_t j = i+1; j < n; ++j) {
				XYZ delta(data[i].x-data[j].x, data[i].y-data[j].y, data[i].z-data[j].z);
				if (shorter(delta, cutoff)) {
					hits += 1;
				}
			}
		}

		Time t2 = now();

		out_log << "Test 6: Same to 4 with buffered size\n"
				"hits " << hits << " and it costed " << t2-t1 << std::endl;
	}

	void TestXYZJob::test7(const std::vector<XYZ>& data, RealNumber cutoff, std::ostream& out_log) {
		Time t1 = now();
		size_t hits = 0;
		RealNumber r2 = cutoff*cutoff;
		for (size_t i = 0; i < data.size(); ++i) {
			for (size_t j = i+1; j < data.size(); ++j) {
				XYZ delta(data[i].x-data[j].x, data[i].y-data[j].y, data[i].z-data[j].z);
				if (shorter2(delta, r2)) {
					hits += 1;
				}
			}
		}

		Time t2 = now();

		out_log << "Test 7: same to 4 with shoter2\n"
				"hits " << hits << " and it costed " << t2-t1 << std::endl;
	}

	void TestXYZJob::test8(const std::vector<XYZ>& data, RealNumber cutoff, std::ostream& out_log) {
		Time t1 = now();
		size_t hits = 0;
		for (size_t i = 0; i < data.size(); ++i) {
			for (size_t j = i+1; j < data.size(); ++j) {
				XYZ delta(data[i], data[j]);
				if (shorter(delta, cutoff)) {
					hits += 1;
				}
			}
		}

		Time t2 = now();

		out_log << "Test 8: Same to 4 with explicit ctor 2\n"
				"hits " << hits << " and it costed " << t2-t1 << std::endl;
	}

	void TestXYZJob::testC1(const std::vector<XYZ>& data, RealNumber cutoff, std::ostream& out_log) {
		Time t1 = now();
		const RealNumber* x = reinterpret_cast<const RealNumber*>(&data[0]);
		size_t n = data.size() * 3;
		size_t hits = 0;
		for (const RealNumber* i = x; i < x+n; i += 3)
		{
			for (const RealNumber* j = i + 3; j < x+n; j += 3) {
				RealNumber delta[3];
				for (size_t k = 0; k < 3; ++k) {
					delta[k] = i[k]-j[k];
				}
				if (fabs(delta[0]) < cutoff && fabs(delta[1]) < cutoff && fabs(delta[2]) < cutoff) {
					RealNumber r2 = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2];
					if (r2 < cutoff*cutoff) {
						hits += 1;
					}
				}
			}
		}
		Time t2 = now();
		out_log << "Test C1: Use C style shorter\n"
				"hits " << hits << " and it costed " << t2-t1 << std::endl;
	}

	void TestXYZJob::testC2(const std::vector<XYZ>& data, RealNumber cutoff, std::ostream& out_log) {
		Time t1 = now();
		const RealNumber* x = reinterpret_cast<const RealNumber*>(&data[0]);
		size_t n = data.size() * 3;
		size_t hits = 0;
		for (const RealNumber* i = x; i < x+n; i += 3)
		{
			for (const RealNumber* j = i + 3; j < x+n; j += 3) {
				RealNumber delta[3];
				for (size_t k = 0; k < 3; ++k) {
					delta[k] = i[k]-j[k];
				}
				if (delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2] < cutoff*cutoff) {
					hits += 1;
				}
			}
		}
		Time t2 = now();
		out_log << "Test C2: Use C style r^2 with pointers\n"
				"hits " << hits << " and it costed " << t2-t1 << std::endl;
	}

	void TestXYZJob::testC3(const std::vector<XYZ>& data, RealNumber cutoff, std::ostream& out_log) {
		Time t1 = now();
		const RealNumber* x = reinterpret_cast<const RealNumber*>(&data[0]);
		size_t n = data.size() * 3;
		size_t hits = 0;
		for (size_t i = 0; i < n; i += 3)
		{
			for (size_t j = i + 3; j < n; j += 3) {
				RealNumber delta[3];
				for (size_t k = 0; k < 3; ++k) {
					delta[k] = x[i+k]-x[j+k];
				}
				if (delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2] < cutoff*cutoff) {
					hits += 1;
				}
			}
		}
		Time t2 = now();
		out_log << "Test C3: Use C style r2, and index\n"
				"hits " << hits << " and it costed " << t2-t1 << std::endl;
	}

	void TestXYZJob::testC4(const std::vector<XYZ>& data, RealNumber cutoff, std::ostream& out_log) {
		Time t1 = now();
		const RealNumber* x = reinterpret_cast<const RealNumber*>(&data[0]);
		size_t n = data.size() * 3;
		size_t hits = 0;
		for (size_t i = 0; i < n; i += 3)
		{
			for (size_t j = i + 3; j < n; j += 3) {
				RealNumber delta[3];
				delta[0] = x[i+0]-x[j+0];
				delta[1] = x[i+1]-x[j+1];
				delta[2] = x[i+2]-x[j+2];
				if (delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2] < cutoff*cutoff) {
					hits += 1;
				}
			}
		}
		Time t2 = now();
		out_log << "Test C4: Use C style r2, and index without inner loop\n"
				"hits " << hits << " and it costed " << t2-t1 << std::endl;
	}

	void TestXYZJob::testC5(const std::vector<XYZ>& data, RealNumber cutoff, std::ostream& out_log) {
		std::vector<RealNumber> x(data.size()), y(data.size()), z(data.size());
		for (size_t i = 0; i < data.size(); ++i) {
			x[i] = data[i].x;
			y[i] = data[i].y;
			z[i] = data[i].z;
		}
		Time t1 = now();
		
		size_t n = data.size();
		size_t hits = 0;
		for (size_t i = 0; i < n; ++i)
		{
			for (size_t j = i + 1; j < n; ++j) {
				RealNumber delta[3] = {x[i]-x[j], y[i]-y[j], z[i]-z[j]};
				if (delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2] < cutoff*cutoff) {
					hits += 1;
				}
			}
		}
		Time t2 = now();
		out_log << "Test C5: split X Y Z\n"
				"hits " << hits << " and it costed " << t2-t1 << std::endl;
	}

	void TestXYZJob::testC6(const std::vector<XYZ>& data, RealNumber cutoff, std::ostream& out_log) {
		std::vector<RealNumber> x(data.size()), y(data.size()), z(data.size());
		for (size_t i = 0; i < data.size(); ++i) {
			x[i] = data[i].x;
			y[i] = data[i].y;
			z[i] = data[i].z;
		}
		RealNumber* px = &x[0];
		RealNumber* py = &y[0];
		RealNumber* pz = &z[0];
		Time t1 = now();
		
		size_t n = data.size();
		size_t hits = 0;
		for (size_t i = 0; i < n; ++i)
		{
			for (size_t j = i + 1; j < n; ++j) {
				RealNumber delta[3] = {px[i]-px[j], py[i]-py[j], pz[i]-pz[j]};
				if (delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2] < cutoff*cutoff) {
					hits += 1;
				}
			}
		}
		Time t2 = now();
		out_log << "Test C6: split X Y Z to C pointers\n"
				"hits " << hits << " and it costed " << t2-t1 << std::endl;
	}


	IMPLEMENT_JOB_CREATOR(TestXYZJob, Test xyz);
	JOB_HELP_FUNCTION(TestXYZJob, output) {
	
	}
}
